The data in linearRegression_crimePunishment.csv contains the murder rate per capita and the rate of automobile crimes per 100,000 individuals (both on the log scale) in the ten US States that have changed their legislation on capital punishment since 1960 (in all cases the states abolished capital punishment). We also include a dummy variable (“law”) that is 1 if the state allows capital punishment in that year, and 0 otherwise. The crime data is from http://www.disastercenter.com.
crime_data <- read.csv('linearRegression_crimePunishment.csv')
crime_data
library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
[30m── [1mAttaching packages[22m ────────────────────────────────────────────────────── tidyverse 1.3.0 ──[39m
[30m[32m✓[30m [34mggplot2[30m 3.3.2 [32m✓[30m [34mpurrr [30m 0.3.4
[32m✓[30m [34mtibble [30m 3.0.3 [32m✓[30m [34mdplyr [30m 1.0.1
[32m✓[30m [34mtidyr [30m 1.1.1 [32m✓[30m [34mstringr[30m 1.4.0
[32m✓[30m [34mreadr [30m 1.3.1 [32m✓[30m [34mforcats[30m 0.5.0[39m
[30m── [1mConflicts[22m ───────────────────────────────────────────────────────── tidyverse_conflicts() ──
[31mx[30m [34mpurrr[30m::[32maccumulate()[30m masks [34mforeach[30m::accumulate()
[31mx[30m [34mdplyr[30m::[32mfilter()[30m masks [34mstats[30m::filter()
[31mx[30m [34mdplyr[30m::[32mlag()[30m masks [34mstats[30m::lag()
[31mx[30m [34mpurrr[30m::[32mwhen()[30m masks [34mforeach[30m::when()[39m
ggplot(data=crime_data) + geom_point(mapping=aes(x=year, y=murder))
ggplot(data=crime_data) + geom_line(mapping=aes(x=year, y=murder))
ggplot(data=crime_data) + geom_smooth(mapping=aes(x=year, y=murder,linetype=stateName))
ggplot(data=crime_data) + geom_smooth(mapping=aes(x=year, y=murder,group=stateName))
ggplot(data=crime_data) + geom_smooth(mapping=aes(x=year, y=murder, color=stateName))
#ggplot(data=filter(crime_data, stateName=='Connecticut')) + geom_line(mapping=aes(x=year, y=murder)) + geom_line(mapping=aes(x=year, y=car))
ggplot(data=filter(crime_data, stateName=='Connecticut'), aes(year)) + geom_line(aes(y=murder, color="murder")) + geom_line(aes(y=car/2, color="car")) + geom_line(aes(y=law, color='law'))
filter(crime_data, stateName=='Connecticut')
df <- filter(crime_data, stateName=='Connecticut') %>%
Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls
3: In readChar(file, size, TRUE) : truncating string with embedded nuls
4: In readChar(file, size, TRUE) : truncating string with embedded nuls
5: In readChar(file, size, TRUE) : truncating string with embedded nuls
6: In readChar(file, size, TRUE) : truncating string with embedded nuls
7: In readChar(file, size, TRUE) : truncating string with embedded nuls
8: In readChar(file, size, TRUE) : truncating string with embedded nuls
select(year, murder, car, law) %>%
gather(key="variable", value = "value", -year)
ggplot(df, aes(x=year, y=value)) +
geom_line(aes(color=variable, linetype=variable)) + ggtitle("Connecticut") + theme(plot.title = element_text(hjust = 0.5))
NA
NA
Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls
3: In readChar(file, size, TRUE) : truncating string with embedded nuls
4: In readChar(file, size, TRUE) : truncating string with embedded nuls
5: In readChar(file, size, TRUE) : truncating string with embedded nuls
6: In readChar(file, size, TRUE) : truncating string with embedded nuls
7: In readChar(file, size, TRUE) : truncating string with embedded nuls
for (lstate in unique(crime_data$stateName)){
df <- filter(crime_data, crime_data$stateName == lstate) %>%
select(year, murder, car, law) %>%
gather(key="variable", value = "value", -year)
print(ggplot(df, aes(x=year, y=value)) +
geom_line(aes(color=variable, linetype=variable)) + ggtitle(lstate) + theme(plot.title = element_text(hjust = 0.5)))
}
A simple model for murder rates is of the form,
murderi,t ∼ N (α + βpenaltyi,t + γcari,t, σ) (18.1)
where we assume that the effect of having the death penalty is given by β, which is assumed to be the same across all states. We include cari,t – a measure of crimes on automobiles, as an independent variable to proxy for the contemporaneous underlying level of crime. Estimate this model and hence determine whether the death penalty acts as a deterrent to murder.
library(rstan)
options(mc.cores = parallel::detectCores())
dataList <- list(N = 540, murder = crime_data$murder, car = crime_data$car,
There were 16 warnings (use warnings() to see them)
law = crime_data$law, state = crime_data$state)
fit <- stan(file = 'Model_Problem_18_1.stan', data = dataList, iter=1000, chains=4, seed=42)
starting worker pid=58687 on localhost:11242 at 17:39:39.666
starting worker pid=58701 on localhost:11242 at 17:39:39.973
starting worker pid=58715 on localhost:11242 at 17:39:40.713
starting worker pid=58734 on localhost:11242 at 17:39:41.209
SAMPLING FOR MODEL 'Model_Problem_18_1' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.000133 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.33 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'Model_Problem_18_1' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 0.000158 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.58 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 1000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'Model_Problem_18_1' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 0.000174 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.74 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 1: Iteration: 100 / 1000 [ 10%] (Warmup)
SAMPLING FOR MODEL 'Model_Problem_18_1' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 0.000169 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.69 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 3: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 2: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 1: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 4: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 3: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 2: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 1: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 4: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 3: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 2: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 1: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 4: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 3: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 2: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 4: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 1: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 1: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 3: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 3: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 2: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 2: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 4: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 4: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 1: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 3: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 1: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 2: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 4: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 3: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 1: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 4: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 2: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 3: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 1: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 4: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 2: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 3: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 1: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 1.16125 seconds (Warm-up)
Chain 1: 1.0339 seconds (Sampling)
Chain 1: 2.19514 seconds (Total)
Chain 1:
Chain 4: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 3: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 1.13128 seconds (Warm-up)
Chain 3: 1.04703 seconds (Sampling)
Chain 3: 2.17831 seconds (Total)
Chain 3:
Chain 2: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 4: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 1.08466 seconds (Warm-up)
Chain 4: 1.16152 seconds (Sampling)
Chain 4: 2.24619 seconds (Total)
Chain 4:
Chain 2: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 1.27196 seconds (Warm-up)
Chain 2: 1.17677 seconds (Sampling)
Chain 2: 2.44873 seconds (Total)
Chain 2:
Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls
3: In readChar(file, size, TRUE) : truncating string with embedded nuls
4: In readChar(file, size, TRUE) : truncating string with embedded nuls
5: In readChar(file, size, TRUE) : truncating string with embedded nuls
6: In readChar(file, size, TRUE) : truncating string with embedded nuls
7: In readChar(file, size, TRUE) : truncating string with embedded nuls
8: In readChar(file, size, TRUE) : truncating string with embedded nuls
print(fit)
Inference for Stan model: Model_Problem_18_1.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
alpha -1.48 0.01 0.22 -1.91 -1.63 -1.48 -1.34 -1.05 784 1
beta 0.24 0.00 0.06 0.12 0.20 0.24 0.28 0.35 1170 1
gamma 0.45 0.00 0.04 0.37 0.43 0.45 0.48 0.53 781 1
sigma 0.62 0.00 0.02 0.59 0.61 0.62 0.64 0.66 1200 1
lp__ -16.05 0.05 1.32 -19.24 -16.77 -15.77 -15.06 -14.38 826 1
Samples were drawn using NUTS(diag_e) at Tue Sep 8 17:39:58 2020.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).